Streszczenie

Obiektem analizy był zbiór danych dotyczących połowu śledzia oceanicznego. Dane zwierały ponad 50 tysięcy wpisów zebranych na przestrzenii około 60 lat. Podczas analizy zbioru zauważono, że od pewnego momentu długość śledzi zaczęła regularnie spadać. Po wstępnym oczyszczeniu danych i usunięciu z nich trywialnych zależności, wykorzystano algorytm uczenia maszynowego eXtreme Gradient Boosting i wykryto, że główną przyczyną tego zjawiska jest utrzymujący się wzrost temperatury przy powierzchni wody.

Użyte biblioteki R

library(dplyr)
library(ggplot2)
library(grid)
library(gridExtra)
library(reshape2)
library(caret)
library(knitr)
library(plotly)

Powtarzalność wyników

set.seed(666)

Zbiór danych

Dane zostały wczytane z pliku CSV. Brakujące wartości oznaczone były znakiem ?, który zamienono na wartość NA.

data <- read.csv("sledzie.csv")
data[data=="?"] <- NA
data <- mutate_all(data, as.numeric)
completeData <- data[complete.cases(data),]

Zbiór danych zawierał 52582 wpisów, każdy z nich opisany był przez 16 atrybutów numerycznych, w tym jedną kolumnę identyfikatora, porządkującą zbiór chronologicznie.

Opis atrybtów

  • length: długość złowionego śledzia [cm],
  • cfin1: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1],
  • cfin2: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2],
  • chel1: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1],
  • chel2: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2],
  • lcop1: dostępność planktonu [zagęszczenie widłonogów gat. 1],
  • lcop2: dostępność planktonu [zagęszczenie widłonogów gat. 2],
  • fbar: natężenie połowów w regionie [ułamek pozostawionego narybku],
  • recr: roczny narybek [liczba śledzi],
  • cumf: łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku],
  • totaln: łączna liczba ryb złowionych w ramach połowu [liczba śledzi],
  • sst: temperatura przy powierzchni wody [°C],
  • sal: poziom zasolenia wody [Knudsen ppt],
  • xmonth: miesiąc połowu [numer miesiąca],
  • nao: oscylacja północnoatlantycka [mb]

Podstawowa charakterystyka

Krótkie podsumowanie najważniejszych statystyk dla każdego z atrybutów (wyłączono atrybut identyfikatora X).

knitr::kable(summary(data[,2:9]))
length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar
Min. :19.0 Min. : 2.00 Min. : 2.00 Min. : 2.00 Min. : 2.00 Min. : 2.00 Min. : 2.00 Min. :0.0680
1st Qu.:24.0 1st Qu.: 2.00 1st Qu.:14.00 1st Qu.:13.00 1st Qu.:15.00 1st Qu.:13.00 1st Qu.:13.00 1st Qu.:0.2270
Median :25.5 Median :15.00 Median :24.00 Median :20.00 Median :27.00 Median :23.00 Median :27.00 Median :0.3320
Mean :25.3 Mean :16.39 Mean :23.64 Mean :23.45 Mean :28.07 Mean :23.85 Mean :27.82 Mean :0.3304
3rd Qu.:26.5 3rd Qu.:28.00 3rd Qu.:34.00 3rd Qu.:36.00 3rd Qu.:41.00 3rd Qu.:35.00 3rd Qu.:43.00 3rd Qu.:0.4560
Max. :32.5 Max. :40.00 Max. :49.00 Max. :49.00 Max. :52.00 Max. :49.00 Max. :52.00 Max. :0.8490
NA NA’s :1581 NA’s :1536 NA’s :1555 NA’s :1556 NA’s :1653 NA’s :1591 NA 
knitr::kable(summary(data[,10:16]))
recr cumf totaln sst sal xmonth nao
Min. : 140515 Min. :0.06833 Min. : 144137 Min. : 2.00 Min. :35.40 Min. : 1.000 Min. :-4.89000
1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068 1st Qu.:14.00 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
Median : 421391 Median :0.23191 Median : 539558 Median :25.00 Median :35.51 Median : 8.000 Median : 0.20000
Mean : 520367 Mean :0.22981 Mean : 514973 Mean :25.45 Mean :35.51 Mean : 7.258 Mean :-0.09236
3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351 3rd Qu.:36.00 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
Max. :1565890 Max. :0.39801 Max. :1015595 Max. :52.00 Max. :35.61 Max. :12.000 Max. : 5.08000
NA NA NA NA’s :1584 NA NA NA

Wstępne przetwarzanie

Brakujące wartości

Pomimo dużej liczby danych w zbiorze, ich zróżnicowanie było bardzo niewielkie. Wynikało to z faktu, że w większości atrybuty opisują wartości mierzone w oknie rocznym (dla danego roku parametry się nie zmieniają). Brakujące wartości w odpowiednich kolumnach zostały zastąpione wartościami z odpowiednich kolumn pełnych rekordów, w których wartość atrybutu recr była taka sama (atrybut ten opisuje roczny narybek, więc dla danego roku jest on stały).

W celu przyśpieszenia wyszukiwania odpowiednich wartości zastępujących, dane zostały pogrupowane wzlędem kolumny recr, a następnie umieszczone w tablicy haszowej. Przedstawiona poniżej funkcja replace_missing zawiera instrukcję warunkową, dla przypadku, gdy dany rekord nie ma atrybutu length - ten przypadek jest zastępowany średnią długością z danego roku i miesiąca. W dostarczonym ziobrze danych, taki przypadek jednak nie wystepował.

replacementData <- group_by(completeData, recr) %>% filter(row_number()==1) %>% ungroup() %>% arrange(recr)
lookupReplacement <- new.env()
for(ridx in 1:nrow(replacementData)){
    key <- as.character(unlist(replacementData[ridx,"recr"]))
    lookupReplacement[[key]] <- replacementData[ridx,]
}

replace_missing <- function(dataRow) {
  recrid <- unlist(dataRow["recr"])
  replacement <- lookupReplacement[[as.character(recrid)]]
  dataRow[which(is.na(dataRow))] <- replacement[which(is.na(dataRow))]
  if(is.na(dataRow["length"])) {
    m <- unlist(dataRow["xmonth"])
    l <- unlist((filter(completeData,recr==recrid,xmonth==m) %>% group_by(recr) %>% summarize(l=mean(length)))["l"])
    dataRow["length"] <- l
  }
  
  
  dataRow
}

withMissing <- data[!complete.cases(data),]
if(nrow(withMissing) > 0){
  for(item in 1:nrow(withMissing)){
    withMissing[item,] <- replace_missing(withMissing[item,])
  }
}

Normalizacja danych

Przed rozpoczęciem analizy statystycznej, dane we wszystkich kolumnach (z wyjątkiem kolumny identyfikatora X) zostały znormalizowane do przedziału od 0.0 do 1.0.

completeData <- arrange(merge(completeData,withMissing,all=TRUE),X) %>% mutate_all(as.numeric)
normalized <- completeData
for(col in colnames(normalized)){
  if(col == "X")next
  normalized <- mutate_(normalized,.dots=setNames(paste0("(",col,"-min(",col,"))/(max(",col,")-min(",col,"))"),col))
}

Korelacja atrybutów

Korelacja pomiędzy atrybutami została sprawdzona za pomocą współczynnika korelacji Pearsona. Wykres przedstawia zależności pomiędzy atrybutami (ciemniejszy kolor, oznacza silniejszą korelację). Dla ułatwienia interpretacji, wykres przedstawia wartości bezwzlędne współczynnika korelacji.

Najsilniesze powiazanie występowało pomiędzy atrybutami:

Długość śledzi

Poniższy animowany wykres przedstawia jak zmieniał się trend długości śledzi z ubiegiem czasu - od początku do końca pomiarów. Łatwo zauważyć, po około 2/5 czasu uwzlędnionego w zbiorze nastąpiło przełamanie trendu i długość śledzi zaczeła maleć.

Wykres został wygenerowany za pomocą kodu w komentarzu, RMarkdown miał problem z użyciem bilbioteki animate i wykres nie generował się poprawnie

Dla łatwiejszej analizy, poniższe wykresy przedstawiają trend dla kilku kluczowych punktów:

library(animation)
show_on<-c(20,40,60,100)
max_x <- max(completeData[,"X"])
x_es <- 1:max(max_x)
x_es <- x_es[x_es%%(floor(length(x_es)/100))==0]
trend_plots <- c()
  for(pseudo_x in x_es){
    percent <- ceiling((pseudo_x*100)/max_x)
    x_label <- paste0("czas, który upłynął: ",percent," %")
    p <- ggplot(filter(completeData, X<=pseudo_x),aes(X,length)) + 
      geom_smooth(method="gam",formula = y ~ s(x, bs = "cs")) + labs(y="trend długości śledzi [cm]",x=x_label) + 
      ylim(25,28) +
      theme( axis.ticks.x = element_blank(), axis.text.x = element_blank())
    trend_plots <- c(trend_plots, p)
    if(any(show_on==percent)){
      print(p)
    }
  }

saveGIF({
  for(idx in seq_along(trend_plots)){
    print(trend_plots[idx])
  }
}, interval = 0.2, ani.width = 550, ani.height = 350)

Rozkład wartości atrybutów

Poniższe histogramy prezentują rozkłady wartości atrybutów w zbiorze. Poza atrybutem length, który charakteryzuje się rozkładem zbliżonym do rozkładu normalnego, atrybuty posiadają bardzo nieregularne rozkłady wartości, z kilkoma wartościami odstającymi. Nie pokazano rozkładu dla atrybutów X (identyfikator) oraz xmonth.

Tak jak wspomniano wcześniej (podpunkt Brakujące Wartości) większość atrybutów zawiera mało unikalnych wartości. Warto zwrócić uwagę na oś Y - liczba elementów, której skala rozciąga się od 0 do 6000. Dla niektórych wartości poszczególnych atrybutów, ich powtarzalność sięga prawie maksymalnej wartości tego przedziału.

Model predykcyjny

Przed rozpoczęciem trenowania modelu, dane zostały znormalizowane do przedziału od 0.0 do 1.0 (wszystkie atrybuty, z wyjątkiem długości - length, który zachował oryginalne wartości). Usunięto również atrybuty o silnej korelacji (na podstawie analizy opisanej wcześniej). Usunięte atrybuty to: chel1, chel2, cumf oraz X (kolumna identyfikatora).

Zbiór został podzielony na treningowy i testujący w stosunku odpowiednio 7:3, przy czym zastosowano losowanie warstowe, zachowujące rozkład wartości atrybutu length. Po wielu eksperymentach, najlepiej ocenianym modelem predykcyjnym okazał się eXtreme Gradient Boosting. Jego parametry zostały dobrane za pomocą metody kroswalidacji.

withoutCorelated <- select(normalized,-chel1,-chel2,-cumf, -X)
herrings <- withoutCorelated
herrings[["length"]] <- completeData[["length"]]
inTraining <- createDataPartition(herrings$length, p=0.7, list=FALSE)
training <- herrings[inTraining,]
testing <- herrings[-inTraining,]
ctrl <- trainControl(method="repeatedcv",number=3,repeats=3)
model <- train(length ~ .,data = training, method="xgbLinear", trControl=ctrl)

Najlepsze parametry oraz miary oceny dla zbioru treningowego:

b<-model$bestTune
kable(b)
nrounds lambda alpha eta
19 50 0.1 0 0.3 
best<-filter(model$results,nrounds==b$nrounds&lambda==b$lambda&alpha==b$alpha&eta==b$eta)
kable(select(best,RMSE,Rsquared))
RMSE Rsquared
1.143814 0.5228055

Na danych testowych model uzyskał następujące miary:

predictions <- testing
predictions[["predicted"]] <- predict(model, newdata=testing)
rmse_testing <- sqrt(mean((predictions$length-predictions$predicted)^2))
ss_res<-sum((predictions$length-predictions$predicted)^2)
ss_tot<-sum((predictions$length-mean(predictions$length))^2)
rsq_testing <- 1 - (ss_res/ss_tot)
kable(data.frame("RMSE"=rmse_testing,"Rsquared"=rsq_testing))
RMSE Rsquared
1.131134 0.5280382 
p <- ggplot(predictions,aes(length,predicted)) + geom_point() + geom_smooth(method="glm") + labs(y="długość przewidziana przez model [cm]",x="rzeczywista długość [cm]") + xlim(19,32.5) + ylim(19,32.5)  + theme_minimal()
ggplotly(p)

Istotność atrybutów

Ranking ważności atrybutów dla zbudowanego modelu regresji przedstawia się następująco:

importance <- varImp(model)
ggplot(importance,aes(Overall)) + geom_bar(stat="identity",fill="#0089B2") + theme_minimal() + labs(y="ważność", x="atrybut",title="Ranking atrybutów")

Najważniejszym atrybutem okazal się sst, oznacza to, że głownym czynnikiem wpływającym na długość śledzia jest temperatura przy powierzchni wody. W początkowym okresie, gdy temperatura spadała, średnia długość śledzi wzrastała. Po osiągnięciu minimum temperatury, zaczęła ona dość gwałtownie rosnąć i od tego samego punktu, długość śledzia zaczęła nieustannie maleć. Na poniższym wykresie można zauwazyć nietrywialną korelację pomiędzy tymi atrybutami:

ggplot(normalized,aes(X,sst)) + geom_smooth(aes(X,sst,colour="temperatura (sst)")) + geom_smooth(mapping=aes(X,length,colour="długość śledzia")) + labs(title="Korelacja sst i length", y="znormalizowana wartość",x="czas") + theme(legend.text =element_text("atrybut")) + scale_colour_manual("",values=c("#0089B2","red")) + theme_minimal() +  theme(legend.position = "bottom")

Pozostałe atrybuty o wiele mniej wpływały na długość śledzi, jednak na niektórych z nich (np. dostępnośc planktonu) można również zauważyć korelację ich trendu względem długości śledzi.

Alternatywny model

Po odkryciu wspomnianych zależności, usunięto 4 najważniejsze atrybuty ze zbioru (sst, recr, totaln, lcop1), aby sprawdzić, jak zachowa się model. Okazało się, że zbudowany model ma bardzo zbliżoną jakość, pomimo wydawałoby się “trudniejszego” zbioru danych. Tym razem, ważność atrybutów jest nieco bardziej zbilansowana - poza najistotniejszym - natężeniu połowów (fbar), atrybuty sal, nao oraz cfin1 są na podobnym poziomie istotności.

herrings2 <- select(withoutCorelated,-lcop1, -sst, -totaln, -recr)
herrings2[["length"]] <- completeData[["length"]]
inTraining2 <- createDataPartition(herrings2$length,p=0.7,list=FALSE)
training2 <- herrings2[inTraining2,]
testing2 <- herrings2[-inTraining2,]
ctrl <- trainControl(method="repeatedcv",number=3,repeats=3)
model2 <- train(length ~ .,data = training2, method="xgbLinear", trControl=ctrl)
predictions2 <- testing2
predictions2[["predicted"]] <- predict(model2, newdata=testing2)
rmse_testing2 <- sqrt(mean((predictions2$length-predictions2$predicted)^2))
ss_res2<-sum((predictions2$length-predictions2$predicted)^2)
ss_tot2<-sum((predictions2$length-mean(predictions2$length))^2)
rsq_testing2 <- 1 - (ss_res2/ss_tot2)
importance2 <- varImp(model2)
ggplot(importance2,aes(Overall)) + geom_bar(stat="identity",fill="#0089B2") + theme_minimal() + labs(x="atrybut", y="ważność",title="Ranking atrybutów modelu alternatywnego")

Wskaźniki dla alternatywnego modelu dla danych treningowych:

RMSE Rsquared
1.146805 0.5211433

oraz testowych:

RMSE Rsquared
1.143814 0.5228055

Podsumowanie

Wejściowy zbiór danych okazał się dość trudny w przetwarzaniu. Dla wielu różnych długości śledzia, wszystkie pozostałe atrybuty posiadały te same wartości, co przyczyniło się do dużych odchyleń w wartościach przewidywanych przez budowane modele predykcyjne.

W trakcie analizy, zgodnie z podanymi informacjami, traktowano X jako atrybut wskazujący na czas, pomimo braku jednostki i jego niespójności ze zmianami wartości atrybutu xmonth. Przy analizie tego typu danych wejściowych, warto poprosić dostawcę danych o przekazanie dodatkowych atrybutów, jeżeli jest taka możliwość.

Wygenerowanie całego raportu od podstaw zajmuje około 5 minut na procesorze i5 @ 1.8GHz